#!/usr/bin/env python
# Copyright 2010 Torbjorn Bjorkman
# This file is part of cif2cell
#
# cif2cell is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# cif2cell is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with cif2cell.  If not, see <http://www.gnu.org/licenses/>.
#
#************************************************************************************
#  Description: A command-line tool to generate the geometrical setup
#               for various electronic structure codes from a CIF format file.
#               This code was published in Comp. Phys. Comm. 182, 1183-1186 (2011).
#  Author:      Torbjorn Bjorkman, torbjorn.bjorkman(at)abo.fi
#  Affiliation: Abo Akademi University
#               Physics/Department of Natural Sciences
#               Porthansgatan 3
#               20500 Turku, Finland
#************************************************************************************
from __future__ import division
import os
import sys
import string
import copy
from math import *
from datetime import datetime
from optparse import OptionParser, OptionGroup
import warnings
import CifFile
import subprocess
from utils import *
from uctools import *
from ESPInterfaces import *
from elementdata import *

# Name and version
programname = "cif2cell"
version = "1.2.10"

# Turn of warnings about deprecated stuff
warnings.simplefilter("ignore",DeprecationWarning)

# All supported output formats
outputprograms = set(['abinit','castep','cfg','coo','cpmd','cp2k','crystal09','elk','emto','exciting','fhi-aims',
                      ## 'fleur','ncol','mcsqs','rspt','siesta','sprkkr','vasp', # mcsqs not ready
                      'fleur','hutsepot','mopac','ncol', 'pwscf', 'quantum-espresso', 'rspt','siesta','sprkkr','vasp',
                      'ase','bmdl','cellgen','cif','kgrn','kfcd','kstr','shape','spacegroup',
                       'spc','xband','xyz','abacus','abacus-lcao'])

# Programs that can deal with alloys
alloyprograms = set(["cfg","emto","kstr","bmdl","shape","kgrn","kfcd","cif","spc","sprkkr","xband"])
vcaprograms = set(["castep","vasp"])
#
setupallprogs = set(["vasp","pwscf","quantum-espresso","rspt","mopac"])
setupallstring = ""
for p in setupallprogs:
    setupallstring += p+", "
setupallstring = setupallstring.rstrip(", ")
#
programlist = sorted(list(outputprograms))
outputprogramstring = ""
for p in programlist:
    outputprogramstring += p+", "
outputprogramstring = outputprogramstring.rstrip(", ")
#print(outputprogramstring)


############################
#      INPUT OPTIONS       #
############################
description = "A program for generating input lattice structures to various electronic structure programs from a CIF (Crystallographic Information Framework) file. This code was published in Comput. Phys. Commun. 182, 1183 (2011). Please cite generously."
usage = "usage: %prog FILE [-p PROGRAM] [other options]"
parser = OptionParser(usage=usage, description=description)
parser.add_option("--version",dest="version",help="Print version number.",action="store_true")
parser.add_option("-v","--verbose",dest="verbose",help="Print a lot of extra information to screen",action="store_true")
parser.add_option("-q","--quiet",dest="quiet",help="Suppress most screen output, such as warnings and overrides --verbose flag, but will print all other requested output.",action="store_true")
# GENERAL OPTIONS
generalopts = OptionGroup(parser, "General options")
generalopts.add_option("-f","--file",dest="file",help="Input CIF file, unless given as first argument to the program.",metavar="FILE")
generalopts.add_option("-p","--program",dest="program",help="The electronic structure code you want to create input file(s) for. Currently supports: "+outputprogramstring+". This keyword is case insensitive.")
generalopts.add_option("-o","--outputfile",dest="outputfile",help="Name of output file (if other than default for you electronic structure code).",metavar="FILE")
generalopts.add_option("-a","--append",dest="append",help="Append the output to given output file rather than overwriting.",action="store_true")
generalopts.add_option("--grammar",dest="grammar",help="Set the CIF grammar to be used when parsing the input file (default is 1.1).")
generalopts.add_option("--which-filename",dest="filenamequery",help="If given together with the --program option, the name of the output file will be printed to screen.",action="store_true")
generalopts.add_option("-b","--block",dest="block",help="Block of data in input file (if there are more than one block in the CIF file).")
# CELL GENERATION OPTIONS
cellgenopts = OptionGroup(parser, "Cell generation options")
cellgenopts.add_option("--no-reduce",dest="noreduce",help="Do not reduce to the primitive cell.",action="store_true")
cellgenopts.add_option("--force",dest="force",help="Attempt to force generation of output file despite problems and/or ambiguities in the input file. There are no guarantees that what you get makes sense, but the program makes an honest attempt. Implies --force-alloy.",action="store_true")
cellgenopts.add_option("--force-alloy",dest="forcealloy",help="Force generation of output file for an alloy compound for an electronic structure code that does not implement any alloy theory (such as CPA).",action="store_true")
cellgenopts.add_option("--vca",dest="vca",help="Set up an alloy using the virtual crystal approximation (VCA). Currently only supported by the CASTEP interface.",action="store_true")
cellgenopts.add_option("--cartesian",dest="cartesian",help="Make the program generate any output in cartesian coordinates.",action="store_true")
cellgenopts.add_option("--coordinate-tolerance",dest="coordtol",help="Parameter for determining when two coordinates are the same (default=0.0002).")
cellgenopts.add_option("--setup-all",dest="setupall",help="Make a more complete setup, not just the geometrical part. This is currently only available for "+setupallstring+".",action="store_true")
cellgenopts.add_option("--kresolution",dest="kresolution",help="The desired resolution in k-space (default=0.2). Used for generating k-space grid options if --setup-all is specified.")
cellgenopts.add_option("--transform-cell",dest="celltransformation",help="Transformation matrix applied to the lattice vectors and the symmetry operations if you, for example, want to realign the cell.",metavar="[[],[],[]]")
cellgenopts.add_option("--cubic-diagonal-z",dest="cubediagz",help="Set up cubic cell with [111] direction along the z-axis.",action="store_true")
cellgenopts.add_option("--rhombohedral-diagonal",dest="rhombdiag",help="Set up rhombohedral cell with threefold axis along pseudocubic [111] direction.",action="store_true")
cellgenopts.add_option("--random-displacements",dest="randomdisp",help="Randomly displace all atoms. Depending on the distribution, the displacement size is either the maximal displacement (for uniform distribution) or the standard deviation (for gaussian distribution) in Angstrom.",metavar="displacementsize")
cellgenopts.add_option("--random-displacements-distribution",dest="randomdistr",help="The distribution used for displacing the atoms.",metavar="uniform/gaussian")
cellgenopts.add_option("--export-cif-labels",dest="exportlabels",help="Export atom labels from the CIF file (currently only supported for castep and RSPt).",action="store_true")
# Supercell options
cellgenopts.add_option("--supercell",dest="supercellmap",help="Three integers separated with commas and enclosed in square brackets that specify the dimensions of a supercell OR three vectors of integers that gives the map to the supercell to be constructed from the primitive cell. If combined with the --no-reduce option the supercell will instead be generated based on the conventional cell.",metavar="[k,l,m]/[[],[],[]]")
cellgenopts.add_option("--supercell-dimensions",dest="supercelldims",help="Three numbers separated with commas and enclosed in square brackets that specify the desired ABSOLUTE dimensions of a supercell (in angstrom) OR three vectors of numbers that gives the desired lattice vectors. The program will automatically generate a supercell, attempting to get as close as possible to the desired dimensions.",metavar="[x,y,z]/[[],[],[]]")
cellgenopts.add_option("--supercell-vacuum",dest="supercellvacuum",help="Three numbers >=0 separated with commas and enclosed in square brackets that specify a number of unit cell units of vacuum to be added along the first, second or third of the generated lattice vectors.",metavar="[k,l,m]")
cellgenopts.add_option("--supercell-translation-vector","--supercell-prevacuum-translation",dest="supercellprevactransvec",help="Three numbers separated with commas and enclosed in square brackets that specify a shift of all atomic positions in the cell prior to vacuum generation (in units of the lattice vectors of the supercell).",metavar="[k,l,m]")
cellgenopts.add_option("--supercell-postvacuum-translation",dest="supercellpostvactransvec",help="Three numbers separated with commas and enclosed in square brackets that specify a final shift of all atomic positions in the final cell (in units of the lattice vectors of the new cell).",metavar="[k,l,m]")
cellgenopts.add_option("--supercell-sort",dest="supercellsort",help="Sort the atom positions by some scheme. Currently available are:     1) By cartesian coordinate - example: xzy will sort first on x then on z then on y.     2) by lattice vector - example: 132 will sort first by lattice vector 1 then by lattice vector 3 and last by lattice vector 2.")
# PRINTING OPTIONS
printopts = OptionGroup(parser, "Printing options")
printopts.add_option("--print-digits",dest="printdigits",help="Number of digits used when printing coordinates etc. to screen (default=8). Useful if you need to tweak the screen output for cutting and pasting into some unsupported program. There is no point in going over 16 because of the floating point accuracy.")
printopts.add_option("--print-atomic-units",dest="printau",help="Output lattice parameters in bohrradii rather than angstrom.",action="store_true")
printopts.add_option("--print-cartesian",dest="printcart",help="Atomic sites printed to screen in cartesian rather than lattice coordinates.",action="store_true")
printopts.add_option("--print-symmetry-operations",dest="printsymops",help="Print symmetry operations of the generated cell.",action="store_true")
printopts.add_option("--print-seitz-matrices",dest="printseitz",help="Print symmetry operations of the generated cell in Seitz matrix form.",action="store_true")
printopts.add_option("--print-charge-state","--print-oxidation-numbers",dest="printcharges",help="Print information about the oxidation state from the CIF file.",action="store_true")
printopts.add_option("--print-reference-bibtex",dest="bibtexref",help="Print citation in BibTeX format and exit.",action="store_true")
# PROGRAM SPECIFIC OPTIONS
progspec = OptionGroup(parser, "Program specific options")
# abinit
progspec.add_option("--abinit-braces",dest="abinitbraces",help="Put curly braces around input values for ABINIT.",action="store_true")
# cellgen
progspec.add_option("--cellgen-map",dest="cellgenmap",help="Nine integers separated with commas and enclosed three and three in square brackets (this is a matrix in Python) that specify the map to a supercell to be output for the RSPt supercell generator 'cellgen'. Overrides --cellgen-supercell-dims.",metavar="[[k,l,m],[n,o,p],[q,r,s]]")
progspec.add_option("--cellgen-supercell-dimensions",dest="cellgensupercelldims",help="Three integers separated with commas and enclosed in square brackets that specify the dimensions of a supercell to be output to the RSPt supercell generator 'cellgen' (the diagonal elements of the 'map').",metavar="[k,l,m]")
progspec.add_option("--cellgen-reference-vector",dest="cellgenrefvec",help="Three reals separated with commas and enclosed in square brackets that specify an optional shift of the origin used by the RSPt supercell generator 'cellgen'.",metavar="[x,y,z]")
# CASTEP
progspec.add_option("--castep-cartesian",dest="castepcartesian",help="Output atom positions in cartesian rather than lattice coordinates.",action="store_true")
progspec.add_option("--castep-atomic-units",dest="castepatomicunits",help="Output to CASTEP in atomic units (bohr radii) rather than angstrom.",action="store_true")
# CPMD
progspec.add_option("--cpmd-cutoff",dest="cpmdcutoff",help="Set the cutoff written to the &SYSTEM block (default=100.0 Ry).")
# emto
progspec.add_option("--emto-hard-sphere-radii",dest="hardsphereradii",help="Set hard spheres in KSTR to something other than the default (=0.67).")
# FHI-AIMS
progspec.add_option("--fhi-aims-cartesian", dest="aimscartesian", help="Store the coordinates for FHI-AIMS in cartesian format.", action="store_true")
# MOPAC
progspec.add_option("--mopac-first-line", dest="mopacfirstline", help="String to be used for the first line (the run commands) of the MOPAC input.",metavar='"string"')
progspec.add_option("--mopac-second-line", dest="mopacsecondline", help="String to be used for the second line (documentation) of the MOPAC input.",metavar='"string"')
progspec.add_option("--mopac-third-line", dest="mopacthirdline", help="String to be used for the third line (documentation) of the MOPAC input.",metavar='"string"')
progspec.add_option("--mopac-freeze-structure", dest="mopacfreeze", help="If set to 'T' then add a 0 after each coordinate (freezing the structure), if set to 'F' then add a 1 (allowing everything to relax).",metavar='T/F')
# PWSCF
progspec.add_option("--pwscf-pseudostring",dest="pwscfpseudostring",help='String to attach to the element name to identify the pseudopotential file (e.g. something like "_HSCV_PBE-1.0.UPF").',metavar="_PSEUDO")
progspec.add_option("--pwscf-atomic-units",dest="pwscfatomicunits",help="Write PWSCF .in file in atomic units (bohr) rather than angstrom.",action="store_true")
progspec.add_option("--pwscf-alat-units",dest="pwscfalatunits",help="Use 'alat' units for the positions in the PWSCF .in file.",action="store_true")
progspec.add_option("--pwscf-cartesian",dest="pwscfcart",help="Write lattice vectors and positions to PWSCF .in file in cartesian coordinates and set the lengths scale to 1.",action="store_true")
progspec.add_option("--pwscf-cartesian-latticevectors",dest="pwscfcartvects",help="Write lattice vectors to PWSCF .in file in cartesian coordinates and set the lengths scale to 1.",action="store_true")
progspec.add_option("--pwscf-cartesian-positions",dest="pwscfcartpos",help="Write lattice positions to PWSCF .in file in cartesian coordinates.",action="store_true")

# ABACUS
progspec.add_option("--deltav",dest="deltav",help='String to attach to lattice volume to modify the original one in cif (e.g. something like "1.1" which means output volume = 1.1 * original volume ).',metavar="1")
progspec.add_option("--energycutoff",dest="energycutoff",help='String to attach to energy cutoff in unit of Ry in output (e.g. something like "80"Ry ).',metavar="80")
progspec.add_option("--orbfilepath",dest="orbfilepath",help='file path to get atmic orbital in lcao method.',metavar="./XXX.dzp")
progspec.add_option("--pseudostringsuffix",dest="pseudostringsuffix",help='String to attach to the element name to identify the pseudopotential file (e.g. something like "_ONCV_PBE-1.0.upf").',metavar="_ONCV_PBE-1.0.upf")
progspec.add_option("--pspdir",dest="pspdir",help='String to attach to the element name to identify the pseudopotential file (e.g. something like "../POT/").',metavar="../POT/")

# RSPt
progspec.add_option("--rspt-new", dest="newsymt", help="Generate a symt.inp file in the new format.", action="store_true")
progspec.add_option("--rspt-spinpol", dest="rsptspinpol", help="Generate new format symt.inp file with spin polarization.", action="store_true")
progspec.add_option("--rspt-relativistic", dest="rsptrelativistic", help="Generate new format symt.inp file with relativistic effects.", action="store_true")
progspec.add_option("--rspt-spinaxis", dest="rsptspinaxis", help="Spin axis for symt.inp (default is [0.0,0.0,0.0].", metavar="[x,y,z]")
progspec.add_option("--rspt-no-spin", dest="rsptnospin", help="Force a nonmagnetic setup in conjunction with --setup-all.", action="store_true")
progspec.add_option("--rspt-mtradii", dest="rsptmtradii", help="Integer that gives the method for setting muffin tin radii.", metavar="N")
progspec.add_option("--rspt-cartesian-latticevectors", dest="rsptcartlatvects", help="Put lattice vectors in atomic units and the lenght scale parameter to 1.", action="store_true")
progspec.add_option("--rspt-pass-wyckoff", dest="rsptpasswyckoff", help="Pass wyckoff labels from CIF file to the symt/rspt.inp file.", action="store_true")
# SPRKKR/xband
progspec.add_option("--sprkkr-minangmom",dest="sprkkrminangmom",help="Enforce minimum onsite angular momentum (=l+1, so that 3 will be d-states).")
# spacegroup
progspec.add_option("--spacegroup-supercell",dest="spacegroupsupercell",help="Three integers separated with commas and enclosed in square brackets that specify the dimensions of a supercell to be output to the elk input generator 'spacegroup'.",metavar="[k,l,m]")
# VASP
progspec.add_option("--vasp-format",dest="vaspformat",help="Format of the generated POSCAR file, either 4 or 5. Default is 4.")
progspec.add_option("--vasp-print-species",dest="vaspprintspcs",help="Print the atomic species to screen in the order they are put in the POSCAR file (useful for scripting).",action="store_true")
progspec.add_option("--vasp-cartesian",dest="vaspcart",help="Write lattice vectors and positions to POSCAR file in cartesian coordinates and set length to 1.",action="store_true")
progspec.add_option("--vasp-cartesian-lattice-vectors",dest="vaspcartvecs",help="Write lattice vectors to POSCAR file in cartesian coordinates and set the length scale to 1.",action="store_true")
progspec.add_option("--vasp-cartesian-positions",dest="vaspcartpos",help="Write atomic positions to POSCAR file in Cartesian rather than Direct coordinates.",action="store_true")
progspec.add_option("--vasp-selective-dynamics",dest="vaspselectivedyn",help="Output POSCAR in selective dynamics format (without any constrained atoms).",action="store_true")
progspec.add_option("--vasp-pseudo-libdr",dest="vasppseudolib",help="Path to the VASP pseudopotential library. Also settable by the VASP_PAWLIB environment variable.")
progspec.add_option("--vasp-pseudo-priority",dest="vasppppriority",help="Set the priority of different pseudopotentials by a list of suffixes. Also available via the VASP_PP_PRIORITY environment variable.", metavar='"_d,_pv,_sv,_h,_s"')
progspec.add_option("--vasp-encutfac",dest="vaspencutfac",help="Factor that multiplies the maximal ENCUT found in the POTCAR file.", metavar="1.5")
# xyz
progspec.add_option("--xyz-atomic-units",dest="xyzatomicunits",help="Output xyz file in atomic units (bohr radii) rather than angstrom.",action="store_true")
#
parser.add_option_group(generalopts)
parser.add_option_group(cellgenopts)
parser.add_option_group(printopts)
parser.add_option_group(progspec)
(options,args) = parser.parse_args()

# Print version number and exit
if options.version:
    print programname+" version "+version
    sys.exit(0)


CellScale  = -1.0
CellVolume = -1.0 
CellNatoms = -1
#global CellVolume
#global CellNatoms
KPTa   = 1
KPTb   = 1
KPTc   = 1
KPTstring = ""

def cif2cellmain():
    global CellScale
    global CellVolume
    global CellNatoms
    global KPTa
    global KPTb
    global KPTc
    global KPTstring

    #############################################################
    # Check that options given are possible
    if options.append and not options.outputfile:
        sys.stderr.write("***Error: option --append requires an output file to be specified.\n")
        sys.exit(1)
    if options.append and (options.program == 'emto' or options.program == 'ncol'):
        sys.stderr.write("***Error: option --append can not be used with "+options.program+".\n")
        sys.exit(1)
    if options.setupall and options.program not in setupallprogs:
        sys.stderr.write("***Error: option --setup-all not supported for "+options.program+".\n")
        sys.exit(1)
    if options.filenamequery and not options.program:
        sys.stderr.write("***Error: option --which-filename requires that --program is given.\n")
        sys.exit(1)
    if options.supercellmap and options.supercelldims:
        sys.stderr.write("***Error: cannot use both --supercell and --supercell-dimensions.")
        sys.exit(1)
    #############################################################
    # INITIAL PARSING OF VARIOUS INPUT DATA
    # Electronic structure program
    if options.program:
        outputprogram = options.program.lower()
        print "    Program:", outputprogram
        if not outputprogram in outputprograms:
            print "Error: Unknown output format: "+outputprogram
            sys.exit(1)
        # Quantum Espresso is just an alias...
        if outputprogram == 'quantum-espresso':
            outputprogram = 'pwscf'
    else:
        outputprogram = None
    
    # recast some input parameters
    if options.noreduce:
        reducetoprim = False
    else:
        reducetoprim = True
    
    try:
        printdigits = int(options.printdigits)
    except:
        printdigits = 8
    
    # Set verbosity level
    if options.verbose and not options.quiet:
        verbose = True
    else:
        verbose = False

    # Various parameters
    # Cartesian?
    if options.cartesian:
        options.castepcartesian = True
        options.printcart = True
        options.vaspcart = True
        options.aimscartesian = True
    # Output reference in specific format?
    if options.bibtexref:
        bibtexref = True
    else:
        bibtexref = False
    # Force generation despite problems?
    if options.force:
        force = True
    else:
        force = False
    # force output for alloys
    if options.forcealloy or force:
        forcealloy = True
    else:
        forcealloy = False
    
    # Make supercell?
    if options.supercellmap or options.supercelldims or options.supercellvacuum or \
           options.supercellprevactransvec or options.supercellpostvactransvec:
        makesupercell = True
    else:
        makesupercell = False
    
    # Initialize element data
    ed = ElementData()
    # Number of positions for printing decimal numbers to screen
    if type(options.printdigits) == type(None):
        decpos = 8 + 3
    else:
        decpos = int(options.printdigits) + 3
    # format string for outputting decimal numbers to screen
    decform = "%"+str(decpos)+"."+str(decpos-4)+"f"
    threedecs = " "+decform+" "+decform+" "+decform
    fourdecs = " "+decform+" "+decform+" "+decform+" "+decform
    # For printing time
    today = datetime.today()
    datestring = str(today.year)+"-"+str(today.month).rjust(2,'0')+"-"+str(today.day).rjust(2,'0')+' '+str(today.hour)+":"+str(today.minute).rjust(2,'0')
    
    # Printing of symmetry operations
    if options.printsymops:
        printsymops = True
    else:
        printsymops = False
    if options.printseitz:
        printseitz = True
    else:
        printseitz = False
    
    if options.printcharges:
        printcharges = True
    else:
        printcharges = False
    
    # complete setup options
    if options.setupall:
        setupall = True
    else:
        setupall = False
    # k-space resolution
    if options.kresolution:
        kresolution=float(options.kresolution)
    else:
        kresolution=0.2
    
    # Cell transformations
    if options.celltransformation or options.cubediagz or options.rhombdiag:
        transformcell = True
    else:
        transformcell = False
        
    #################################################################
    # Open and read CIF file
    cif_file = None
    if len(args) > 0:
        # input CIF file as argument
        cif_file = args[0]
    if options.file:
        # input CIF file as option (overrides argument)
        cif_file = options.file
    if cif_file:
        if not os.path.exists(cif_file):
            sys.stderr.write("***Error: The file "+cif_file+" could not be found.\n")
            sys.exit(2)
        cif_file_name = cif_file.split("/")[-1]
        # Set CIF grammar
        if options.grammar:
            cif_grammar = options.grammar
        else:
            cif_grammar = '1.1'
        # Skip validation for now... it causes too much trouble.
        ## cdic = CifFile.CifDic("cif_core.dic",grammar='1.1')     
        ## val_results = CifFile.validate(cif_file,dic=cdic)
        ## print validate_report(val_results)
        ## val_report = CifFile.ValidationResult(val_results)
        try:
            cf = CifFile.ReadCif(cif_file,grammar=cif_grammar)
        except Exception, e:
            # test if data_ statement in the beginning is missing
            try:
                f = open(cif_file,'r')
                lines = f.readlines()
                f.close()
                tmpname = cif_file.replace('.cif','_tmp.cif') 
                f = open(tmpname,'w')
                f.write("data_default\n")
                for line in lines:
                    f.write(line)
                f.close()
                cf = CifFile.ReadCif(tmpname,grammar=cif_grammar)
                wrongfile = cif_file.replace('.cif','_wrong.cif')
                sys.stderr.write("***Warning: The cif file is missing a data statement.")
                sys.stderr.write(" The file has been renamed '"+wrongfile+"' and replaced by a")
                sys.stderr.write(" corrected file.\n")
                os.rename(cif_file,wrongfile)
                os.rename(tmpname,cif_file)
            except:
                try:
                    os.remove(tmpname)
                except:
                    pass
                sys.stderr.write("***Error: could not read "+cif_file+".\n")
                sys.stderr.write("Something may be wrong with the CIF file, you can check it with ")
                sys.stderr.write("the free IUCr CIF valitation tool at http://http://checkcif.iucr.org/\n")
                sys.stderr.write(e.value+"\n")
                sys.exit(2)
    else:
        sys.stderr.write("***Error: No input CIF file given\n")
        sys.exit(2)
    
    # Make supercell?
    if makesupercell:
        if options.supercellmap:
            supercellmap = safe_matheval(options.supercellmap)
        else:
            supercellmap = [1,1,1]
        if options.supercelldims:
            t = safe_matheval(options.supercelldims)
            try:
                supercelldims = [[float(t[0]), 0.0, 0.0],
                                 [0.0, float(t[1]), 0.0],
                                 [0.0, 0.0, float(t[2])]]
            except:
                supercelldims = t
        else:
            supercelldims = None
        if options.supercellvacuum:
            supercellvacuum = safe_matheval(options.supercellvacuum)
        else:
            supercellvacuum = [0,0,0]
        if options.supercellprevactransvec:
            supercellprevactransvec = safe_matheval(options.supercellprevactransvec)
        else:
            supercellprevactransvec = [0,0,0]
        if options.supercellpostvactransvec:
            supercellpostvactransvec = safe_matheval(options.supercellpostvactransvec)
        else:
            supercellpostvactransvec = [0,0,0]
        # 
        if options.supercellsort:
            supercellsort = options.supercellsort.lower()
        else:
            supercellsort = ""
    
    ##############################################
    # Get blocks
    cfkeys = cf.keys()
    if options.block:
        cb = cf.get(options.block)
        if type(cb) == type(None):
            sys.stderr.write("***Error: No block "+options.block+" in "+cif_file+".\n")
            sys.exit(2)
    else:
        cb = cf.get(cfkeys[0])
    # Get reference data
    ref = ReferenceData()
    ref.getFromCIF(cb)
    if bibtexref:
        print ref.bibtexref()
        sys.exit(0)
    # Get cell data
    cd = CellData()
    # Suppress warnings if requested.
    cd.quiet = options.quiet
    # Force generation despite problems?
    cd.force = force
    if options.coordtol:
        cd.coordepsilon = float(options.coordtol)
    try:
        cd.getFromCIF(cb)
    except PositionError, e:
        sys.stderr.write("***Error: cell setup: "+e.value+"\n")
        sys.exit(2)
    except CellError, e:
        sys.stderr.write("***Error: cell setup: "+e.value+"\n")
        sys.exit(2)
    except SymmetryError, e:
        sys.stderr.write("***Error: cell setup: "+e.value+"\n")
        sys.exit(2)
    
    
    ##############################################
    # Generate cell
    try:
        if reducetoprim:
            cd.primitive()
        else:
            cd.conventional()
    except SymmetryError, e:
        sys.stderr.write("***Error: cell setup: "+e.value+"\n")
        sys.exit(2)
    except CellError, e:
        sys.stderr.write("***Error: cell setup: "+e.value+"\n")
        sys.exit(2)
    
    # Test if generated cell agrees with given chemical formula.
    # Too difficult for alloys.
    if len(ref.ChemicalComposition) > 0 and not cd.alloy:
        if ref.ChemicalComposition != cd.ChemicalComposition:
            if force:
                sys.stderr.write("***Warning: Chemical composition of the generated cell differs from that given\n"+\
                    "            by _chemical_formula_sum.\n")
            else:
                sys.stderr.write("***Error: Chemical composition of the generated cell differs from that given\n"+\
                    "          by _chemical_formula_sum. Use --force to generate a cell anyway.\n")
                sys.exit(2)
    
    inputcell = copy.copy(cd)
    # Randomly displace atoms if requested. This erases all symmetry operations.
    if options.randomdisp and not makesupercell:
        if options.randomdistr:
            distr = options.randomdistr
        else:
            distr = "uniform"
        try:
            cd.randomDisplacements(float(options.randomdisp),distribution=distr)
        except SetupError, e:
            sys.stderr.write("***Error: random displacements: "+e.value+"\n")
            sys.exit(3)
        # Reset space group operations
        cd.HallSymbol = "P 1"
        cd.spacegroupnr = 1
        cd.HMSymbol = "P1"
        cd.symops = set([SymmetryOperation(['x','y','z'])])
    
    # Print cell
    if verbose or not options.program and not options.quiet:
        print string.upper(programname)+" "+version
        print datestring
        # Print compound
        compoundstring = "Output for "
        if ref.cpd == "" and ref.compound == "":
            compoundstring += "unknown compound"
        if ref.cpd != "":
            compoundstring += ref.cpd
        if ref.compound != "":
            compoundstring += " ("+ref.compound+")"
        print compoundstring
        # Print database
        print ref.databasestring
        if cd.alloy and forcealloy and options.program:
            print "\nEnforcing generation of file(s) for "+outputprogram+" for an alloy."
        print "\n BIBLIOGRAPHIC INFORMATION"
        refstrings = ref.referencestring().split()
        tmpstring = ""
        i = 0
        while i < len(refstrings):
            if len(tmpstring+refstrings[i]+" ") < 70:
                tmpstring += refstrings[i]+" "
                i += 1
            else:
                print tmpstring
                tmpstring = ""
        if tmpstring != "":
            print tmpstring
        print "\n INPUT CELL INFORMATION"
        print "Symmetry information:"
        if inputcell.HallSymbol != "":
            print inputcell.crystal_system()[0].upper()+inputcell.crystal_system()[1:]+" crystal system."
            print "Space group number     : ".rjust(2)+str(inputcell.spacegroupnr)
            print "Hall symbol            : "+inputcell.HallSymbol
            print "Hermann-Mauguin symbol : "+inputcell.HMSymbol
        else:
            print "No space group information found."
        # only print these if verbose
        if verbose:
            print "Symmetry equivalent sites:"
            symops = list(inputcell.symops)
            symops.sort()
            for i in range(len(symops)):
                print "%4i  %8s, %8s, %8s" % (i+1, symops[i].eqsite[0], symops[i].eqsite[1], symops[i].eqsite[2])
        print "\nLattice parameters:"
        tmpstring = ""
        for i in ["a", "b", "c"]:
            tmpstring += i.rjust(decpos)+" "
        print tmpstring
        formatstring = ""
        if options.printau:
            ## aprint = inputcell.ainit*angtobohr
            ## bprint = inputcell.binit*angtobohr
            ## cprint = inputcell.cinit*angtobohr
            aprint = inputcell.a*angtobohr
            bprint = inputcell.b*angtobohr
            cprint = inputcell.c*angtobohr
        else:
            ## aprint = inputcell.ainit
            ## bprint = inputcell.binit
            ## cprint = inputcell.cinit
            aprint = inputcell.a
            bprint = inputcell.b
            cprint = inputcell.c
        for i in range(3):
            formatstring = formatstring+decform+" "
        print formatstring % (aprint, bprint, cprint)
        tmpstring = ""
        for i in ["alpha", "beta", "gamma"]:
            tmpstring += i.rjust(decpos)+" "
        print tmpstring
        print formatstring % (inputcell.alpha, inputcell.beta, inputcell.gamma)
        ## print formatstring % (inputcell.alphainit, inputcell.betainit, inputcell.gammainit)
        # Pretty printing in columns that need to have variable width
        # w1 = width of the atomic species column
        # w2 = width of a decimal column
        # w3 = width of the occupancy column
        # w4 = width of the charge state column
        if inputcell.alloy:
            w1 = 0
            w3 = 0
            w4 = 0
            # Find atom and occupation column widths
            for a in inputcell.atomdata:
                for b in a:
                    tmpstring1 = ""
                    tmpstring2 = ""
                    tmpstring3 = ""
                    for k,v in b.species.iteritems():
                        tmpstring1 += k+"/"
                        tmpstring2 += str(v).rstrip("0.")+"/"
                        # charge output
                        for k2,v2 in inputcell.chargedict.iteritems():
                            if k2.strip(string.punctuation+string.digits) == k:
                                tmpstring3 += str(v2)+"/"
                    tmpstring1 = tmpstring1.rstrip("/")
                    tmpstring2 = tmpstring2.rstrip("/")
                    tmpstring3 = tmpstring3.rstrip("/")
                    w1 = max(w1,len(tmpstring1))
                    w3 = max(w3,len(tmpstring2))
                    w4 = max(w4,len(tmpstring3))
            # small aesthetic adjustment 
            w1 = w1 + 1
            w3 = w3 + 2
            w4 = max(w4 + 2, 8)
        else:
            w1 = 5
            w2 = decpos
            w3 = 0
            # width of charge column
            if printcharges:
                w4 = 7
            else:
                w4 = 0
        # Now for the output...
        tmpstring = "Representative sites :"
        print tmpstring
        siteheader = "Atom".ljust(w1)+" "
        if options.printcart:
            transmtx = []
            for i in range(3):
                transmtx.append([])
                for j in range(3):
                    transmtx[i].append(inputcell.latticevectors[i][j]*inputcell.lengthscale)
                i += 1
            for i in ["x","y","z"]:
                siteheader += i.rjust(decpos)+" "
        else:
            transmtx = [[1, 0, 0],
                        [0, 1, 0],
                        [0, 0, 1]]
            for i in ["a1","a2","a3"]:
                siteheader += i.rjust(decpos)+" "
        if inputcell.alloy:
            if w3 > 13:
                siteheader += "occupancies".rjust(w3)
            else:
                siteheader += "occ.".rjust(w3)
        if printcharges:
            siteheader += " "+"charge".rjust(w4)
        print siteheader
        # Representative sites
        for i in range(len(inputcell.ineqsites)):
            tmpstring = ""
            occstring = ""
            chargestring = ""
            for k,v in inputcell.occupations[i].iteritems():
                tmpstring += k+"/"
                occstring += str(v)+"/"
                # charge output
                for k2,v2 in inputcell.chargedict.iteritems():
                    if k2.strip(string.punctuation+string.digits) == k:
                        chargestring += str(v2)+"/"
            tmpstring = tmpstring.rstrip("/")
            occstring = occstring.rstrip("/")
            chargestring = chargestring.rstrip("/")
            v = [t for t in inputcell.ineqsites[i]]
            tmpstring = tmpstring.ljust(w1) + threedecs % (v[0],v[1],v[2])
            if inputcell.alloy:
                tmpstring += " "+occstring.rjust(w3)
            if printcharges:
                tmpstring += " "+chargestring.rjust(w4)
            print tmpstring
    
        # Output cell
        print "\n OUTPUT CELL INFORMATION"
        print "Symmetry information:"
        if cd.HallSymbol != "":
            print cd.crystal_system()[0].upper()+cd.crystal_system()[1:]+" crystal system."
            print "Space group number     : ".rjust(2)+str(cd.spacegroupnr)
            print "Hall symbol            : "+cd.HallSymbol
            print "Hermann-Mauguin symbol : "+cd.HMSymbol
        else:
            print "No space group information found."
        # only print these if verbose
        if verbose:
            print "Symmetry equivalent sites:"
            symops = list(cd.symops)
            symops.sort()
            for i in range(len(symops)):
                print "%4i  %8s, %8s, %8s" % (i+1, symops[i].eqsite[0], symops[i].eqsite[1], symops[i].eqsite[2])
    
        print ""
        cd.printCell(printcart=options.printcart, printdigits=printdigits, printcharges=options.printcharges)
        # Print volume and density
        if options.printau:
            volume = cd.volume()*(cd.lengthscale*angtobohr)**3
            volstring = "(a.u.)"
        else:
            volume = cd.volume()*cd.lengthscale**3
            volstring = "A"
        print "\nUnit cell volume  : "+decform%volume+" "+volstring+"^3"
        try:
            weight = 0.0
            for a in cd.atomdata:
                for b in a:
                    for k,v in b.species.iteritems():
                        weight += ed.elementweight[k]*v
            density = weight/volume
            print "Unit cell density : "+decform%density+" u/"+volstring+"^3 = "+decform%(density*uperautogpercm)+" g/cm^3"
        except:
            if not options.quiet:
                sys.stderr.write("***Warning: Error printing unit cell density.\n")
    
    ##############################################
    # Rotate cell
    if transformcell:
        # pre-defined transformations
        if options.cubediagz:
            # Put [111] direction along z axis.
            celltransformation = LatticeMatrix([[0.577350269189626,-1.000000000000000,0.816496580927725],
                                                [0.577350269189626, 1.000000000000000,0.816496580927725],
                                                [-1.154700538379250,0.000000000000000,0.816496580927725]])
            if cd.crystal_system() != 'cubic':
                sys.stderr.write("***Error: Only cubic structures are properly aligned to the z axis by --cubic-diagonal-z.\n")
                if not force:
                    sys.stderr.write("          Use --force to go ahead anyway.\n")
                    sys.exit(1)
                
        if options.rhombdiag:
            # Put z direction along pseudocubic [111]. Kind of the opposite of the cubediagz.
            t = 1/cd.latticevectors[0].length()   # Normalization to 1
            celltransformation = LatticeMatrix([[0.816496580927726*t,-0.408248290463863*t,-0.408248290463863*t],
                                                [0.000000000000000, 0.707106781186547*t,-0.707106781186547*t],
                                                [0.577350269189626*t, 0.577350269189626*t, 0.577350269189626*t]])
            if not cd.rhombohedral:
                sys.stderr.write("***Error: Only rhombohedral cells are properly aligned to the\n")
                sys.stderr.write("          pseudocubic (111) axis by --rhombohedral-diagonal.\n")
                if not force:
                    sys.stderr.write("          Use --force to go ahead anyway.\n")
                    sys.exit(1)
        # explicitly giving the transformation overrides any other transformation
        if options.celltransformation:
            celltransformation = safe_matheval(options.celltransformation)
        try:
            cd.transformCell(celltransformation)
        except CellError, e:
            sys.stderr.write("***Error: Cell transformation: "+e.value+"\n")
            sys.exit(2)
        if verbose or not options.program and not options.quiet:
            print "\n TRANSFORMED CELL"
            cd.printCell(printcart=options.printcart, printdigits=printdigits, printcharges=options.printcharges)
        
    ##############################################
    # Generate supercell
    if makesupercell:
        if supercelldims != None:
            # Determine a suitable map to get the desired supercell dimensions.
            t1 = []
            for i in range(3):
                t1.append([])
                for j in range(3):
                    t1[i].append(cd.latticevectors[i][j]*cd.lengthscale)
            t2 = minv3(t1)
            t2 = mmmult3(supercelldims,t2)
            supercellmap = []
            for i in range(3):
                supercellmap.append([])
                for j in range(3):
                    supercellmap[i].append(int(round(t2[i][j])))
        try:
            cd.getSuperCell(supercellmap,supercellvacuum,supercellprevactransvec,postvactransvec=supercellpostvactransvec,sort=supercellsort)
        except CellError, e:
            sys.stderr.write("***Error: Supercell setup: "+e.value+"\n")
            sys.exit(2)
        # Randomly displace atoms if requested. This erases all symmetry operations.
        if options.randomdisp:
            if options.randomdistr:
                distr = options.randomdistr
            else:
                distr = "uniform"
            try:
                cd.randomDisplacements(float(options.randomdisp),distribution=distr)
            except SetupError,e:
                sys.stderr.write("***Error: random displacements: "+e.value+"\n")
                sys.exit(3)
            cd.symops = set([SymmetryOperation(['x','y','z'])])
        
        # Print supercell
        if verbose or not options.program and not options.quiet:
            print "\n SUPERCELL INFORMATION"
            cd.printCell(printcart=options.printcart, printdigits=printdigits, printcharges=options.printcharges)
            
    if printsymops or printseitz or verbose:
        # Print symmetry operations. Need to make list of it to control order.
        symoplist = sorted(list(cd.symops))
        if printsymops or verbose:
            print "\nSymmetry operations : "
            print "  3x3 rotation matrix +"
            print "  3x1 translation vector"
            i = 1
            for op in symoplist:
                print "Operation "+str(i)
                for v in op.rotation:
                    print threedecs%(v[0],v[1],v[2])
                print threedecs%(op.translation[0],op.translation[1],op.translation[2])
                i += 1
        if printseitz:
            print "\nSymmetry operations :"
            print "  In Seitz matrix form"
            i = 1
            for op in symoplist:
                print "Operation "+str(i)
                tmpstring = ""
                for j in range(3):
                    tmpstring += fourdecs%(op.rotation[j][0],op.rotation[j][1],op.rotation[j][2],op.translation[j])+"\n"
                tmpstring += fourdecs%(0,0,0,1)
                print tmpstring
                i += 1    
    
    # Remind that the result may be junk when using --force
    if force:
        sys.stderr.write("\n***Warning: You invoked the --force flag, presumably to bypass some error message.\n")
        sys.stderr.write("            Carefully check the results, which may be rubbish, nonsense or both!\n")
    
    ##############################################
    # Sort sites so that the ones occupied by the heaviest elements come first,
    # if the Python version supports this form of the max function and there is
    # nothing wrong with the site data.
    if not (makesupercell and supercellsort):
        try:
            cd.atomdata.sort(key = lambda a: ed.elementnr[max(a[0].species, key = a[0].species.get)], reverse=True)
        except:
            pass
    
    ############################################################################################
    # Output file mode (overwrite or append?)
    if options.append:
        outmode = "a"
    else:
        outmode = "w"
    # Output file. ot parsed until this point, since the default file names for some of
    # the codes contain the names of space group and compound.
    if options.outputfile:
        outputfile = options.outputfile
    else:
        # Default output filenames for different programs
        if outputprogram == "vasp":
            outputfile = "POSCAR"
        elif outputprogram == "rspt":
            if setupall:
                outputfile = "rspt.inp"
            else:
                outputfile = "symt.inp"
        elif outputprogram == "cellgen":
            outputfile = "cellgen.inp"
        elif outputprogram == "elk":
            outputfile = "GEOMETRY.OUT"
        elif outputprogram == "abacus":
            outputfile = "INPUT"
        elif outputprogram == "abacus-lcao":
            outputfile = "INPUT"
        elif outputprogram == "exciting":
            outputfile = "input.xml"
        elif outputprogram == "spacegroup":
            outputfile = "spacegroup.in"
        elif outputprogram == "cif":
            outputfile = cif_file_name.replace(".cif","")+"_allatoms.cif"
        elif outputprogram == "ase":
            outputfile = "positions.py"
        else:
            # A bunch of programs get default output constructed as:
            # 1. chemical abbreviation (i.e. something like "H2SO4" or CeRhIn5)
            # 2. if this is too long, use the original cif filename
            outputfile = ref.cpd.replace(" ", "").replace("(","").replace(")","")
            # If the filename seems too long or strange, replace by the name of the cif file
            if len(outputfile.strip(string.punctuation)) == 0:
                outputfile = cif_file_name.replace(".cif","")
            if len(outputfile) > 24:
                if len(cif_file_name) < len(outputfile):
                    outputfile = cif_file_name.replace(".cif","")
                else:
                    outputfile = outputfile[0:9]
            # Append file endings etc.
            if outputprogram == "abinit":
                outputfile = outputfile+".in"
            elif outputprogram == "castep":
                outputfile = outputfile+".cell"
            elif outputprogram == "cfg":
                outputfile = outputfile+".cfg"
            elif outputprogram == "coo":
                outputfile = outputfile+".coo"
            elif outputprogram == "cp2k":
                outputfile = outputfile+".inp"
            elif outputprogram == "cpmd":
                outputfile = outputfile+".inp"
            elif outputprogram == "crystal09":
                # This is the naming convention from a large bunch of test cases, no idea why
                outputfile = outputfile+".d12"
            elif outputprogram == "fhi-aims":
                outputfile = "geometry.in"
            elif outputprogram == "fleur":
                outputfile = "inp_"+outputfile
            elif outputprogram == "hutsepot":
                outputfile = outputfile+".str"
            elif outputprogram == "mopac":
                outputfile = outputfile+".mop"
            elif outputprogram == "pwscf":
                outputfile = outputfile+".in"
            elif outputprogram == "siesta":
                outputfile = outputfile+".fdf"
            elif outputprogram == "sprkkr" or outputprogram == "xband":
                outputfile = outputfile+".sys"
            elif outputprogram == "spc":
                outputfile = outputfile+".dat"
            elif outputprogram == "xyz":
                outputfile = outputfile+".xyz"
            else:
                pass
    
    # Print output filename to screen
    if outputprogram !="":
        if (verbose or options.filenamequery) and outputfile != "":
            print "Data will be written to the file "+outputfile
    
    ################################################################################################
    # stuff that should be printed irrespective of the verbose flag
    if cd.alloy and forcealloy and options.program and not verbose:
        tmpstring = "Enforcing file generation for alloy. "
        if outputprogram == 'bstr' or outputprogram == 'vasp' or outputprogram == 'cpmd':
            print tmpstring
        else:
            tmpstring += "Warning! The file(s) will be incomplete!"
            print tmpstring
        
    ################################################################################################
    # Stop here if no specific output was requested
    if not options.program:
        sys.exit(0)
    
    # Don't generate output for alloys (for most programs)
    if cd.alloy and not forcealloy and not (outputprogram in alloyprograms or outputprogram in vcaprograms):
        print "Error: This system is an alloy, but "+codename[outputprogram]+" has no way of dealing with alloys.\n       Run again with --force-alloy (or --force) if you want to generate an (incomplete) output file anyway."
        sys.exit(17)
    # Deal with VCA
    if cd.alloy and outputprogram in vcaprograms:
        if not options.vca and not forcealloy:
            print "Error: This system is an alloy. "+codename[outputprogram]+" can deal with some alloys using the virtual crystal approximation (VCA).\n       Run again with the flag --vca if you want to produce a VCA setup."
            sys.exit(17)
    vcawarning1 = False
    vcawarning2 = False
    if options.vca:
        # Issue warning for precarious VCA setups
        groups = []
        for a in cd.atomdata:
            if len(a[0].species) > 1:
                t = [ed.elementgroup[sp] for sp,conc in a[0].species.iteritems()]
                groups.append((min(t),max(t)))
                if len(a[0].species) > 2:
                    vcawarning1 = True
        for g in groups:
            if g[1] - g[0] > 1:
                vcawarning2 = True
        if vcawarning1 and vcawarning2:
            sys.stderr.write("Warning: You are setting up a VCA calculation for an alloy with more than two components\n         and not all alloy sites are occupied by species from neighbouring groups in the periodic\n         table. Make doubly sure that you know what you are doing!\n")
        elif vcawarning1:
            sys.stderr.write("Warning: You are setting up a VCA calculation for an alloy with more than two components.\n         Make sure that you know what you are doing!\n")
        elif vcawarning2:
            sys.stderr.write("Warning: You are setting up a VCA calculation but not all alloy sites are occupied by species\n         from neighbouring groups in the periodic table. Make sure that you know what you are doing!\n")
            
    # Function for printing a standard docstring
    def StandardDocstring():
        cif2cellstring = ' T. Bjorkman, Comp. Phys. Commun. 182, 1183-1186 (2011). Please cite generously.'
        stringlen = max(len(ref.referencestring()),len(ref.cpd+"   ("+ref.compound+")"),len(cif2cellstring))
        docstring = ""
        tmpstring = ""
        tmpstring = tmpstring.ljust(stringlen+4,'*')+'\n'
        docstring += tmpstring
        tmpstring2 = 'Generated by '+programname+' '+version+' '+datestring
        tmpstring2 = '* '+tmpstring2.center(stringlen)+' *\n'
        tmpstring3 = '* '+cif2cellstring.center(stringlen)+' *\n'
        tmpstring4 = ''
        tmpstring4 = '* '+tmpstring4.center(stringlen)+' *\n'
        docstring += tmpstring2+tmpstring3+tmpstring4
        if ref.database != "":
            tmpstring2 = 'Data obtained from '+ref.databaseabbr[ref.database]
            if ref.databasecode != "":
                tmpstring2 += ". Reference number : "+ref.databasecode
            tmpstring2 = '* '+tmpstring2.center(stringlen)+' *\n'
            docstring += tmpstring2
        tmpstring2 = ref.cpd+"   ("+ref.compound+")"
        tmpstring2 = '* '+tmpstring2.center(stringlen)+' *\n'
        docstring += tmpstring2
        docstring += '* '+ref.referencestring().center(stringlen)+' *\n'
        docstring += tmpstring
        return docstring
    
    ################################################################################################
    # Output cell to new CIF file
    if outputprogram == 'cif':
        f = open(outputfile,'w')
        cf = CifFile.CifFile()
        cb = CifFile.CifBlock()
        if makesupercell or (reducetoprim and cd.spacegroupsetting != 'P') or options.randomdisp:
            a = Vector(cd.latticevectors[0].scalmult(cd.lengthscale))
            b = Vector(cd.latticevectors[1].scalmult(cd.lengthscale))
            c = Vector(cd.latticevectors[2].scalmult(cd.lengthscale))
            cb['_cell_length_a'] = a.length()
            cb['_cell_length_b'] = b.length()
            cb['_cell_length_c'] = c.length()
            cb['_cell_angle_alpha'] = acos(b.dot(c)/(b.length()*c.length()))*180/pi
            cb['_cell_angle_beta'] = acos(a.dot(c)/(a.length()*c.length()))*180/pi
            cb['_cell_angle_gamma'] = acos(b.dot(a)/(a.length()*b.length()))*180/pi
            # Supercell may have broken symmetry, so just put P1
            cb['_space_group_IT_number'] = 1
            cb['_space_group_name_H-M_alt'] = "P1"
            cb['_space_group_name_Hall'] = "P 1"
        else:
            # Else pass on original cell parameters and symmetry information
            cb['_cell_length_a'] = cd.a
            cb['_cell_length_b'] = cd.b
            cb['_cell_length_c'] = cd.c
            cb['_cell_angle_alpha'] = cd.alpha
            cb['_cell_angle_beta']  = cd.beta
            cb['_cell_angle_gamma'] = cd.gamma
            cb['_space_group_IT_number'] = cd.spacegroupnr
            cb['_space_group_name_H-M_alt'] = cd.HMSymbol
            cb['_space_group_name_Hall'] = cd.HallSymbol
        # Positions
        labels = []
        symbols = []
        symmult = []
        fractx = []
        fracty = []
        fractz = []
        occup = []
        i = 0
        for a in cd.atomdata:
            if not makesupercell:
                i += 1
            for b in a:
                if makesupercell:
                    i += 1
                for k in b.species:
                    labels.append(k+str(i))
                    symbols.append(k)
                    if makesupercell:
                        symmult.append("1")
                    else:
                        symmult.append(str(len(a)))
                    fractx.append(("%19.16f"%b.position[0]).strip(" "))
                    fracty.append(("%19.16f"%b.position[1]).strip(" "))
                    fractz.append(("%19.16f"%b.position[2]).strip(" "))
                    occup.append(str(b.species[k]))
        #
        cb.AddCifItem(([['_atom_site_label',
                         '_atom_site_type_symbol',
                         '_atom_site_symmetry_multiplicity',
                         '_atom_site_fract_x',
                         '_atom_site_fract_y',
                         '_atom_site_fract_z',
                         '_atom_site_occupancy']],
                       [[labels,symbols,symmult,fractx,fracty,fractz,occup]]))
        #
        cf['1-cif2cell'] = cb
        f.write(str(cf))
        f.close()
    
    ################################################################################################
    # Output for ABINIT
    if outputprogram == 'abinit':
        docstring = StandardDocstring()
        abinitinput = ABINITFile(cd, docstring)
        if options.abinitbraces:
            abinitinput.printbraces = True
        f = open(outputfile, outmode)
        f.write(str(abinitinput))
        f.close()
        
    ################################################################################################
    # Output for ASE file (python script)
    if outputprogram == 'ase':
        # The second comment line with info from the CIF file 
        docstring = StandardDocstring()
        # Initialize the ASEFile structure
        sysfile = ASEFile(cd, docstring)
        f = open(outputfile, outmode)
        f.write(str(sysfile))
        f.close()
        sys.exit(0)
    
    ################################################################################################
    # Output for CASTEP
    if outputprogram == 'castep':
        docstring = StandardDocstring()
        castepinput = CASTEPFile(cd, docstring)
        if options.castepatomicunits:
            castepinput.unit = "bohr"
        if options.castepcartesian:
            castepinput.cartesian = True
        if options.vca:
            castepinput.vca = True
        if options.exportlabels:
            castepinput.printlabels = True
        f = open(outputfile, outmode)
        f.write(str(castepinput))
        f.close()
        
    ################################################################################################
    # Output to cfg file
    if outputprogram == 'cfg':
        docstring = StandardDocstring()
        # Initialize the CFGFile structure
        sysfile = CFGFile(cd, docstring)
        f = open(outputfile, outmode)
        f.write(str(sysfile))
        f.close()
        sys.exit(0)
    
    ################################################################################################
    # Output to coo file
    if outputprogram == 'coo':
        # !!!TODO: Think of something better for docstring !!!
        docstring = "Generated by cif2cell "+version
        # Initialize the COOFile structure
        sysfile = COOFile(cd, docstring)
        f = open(outputfile, outmode)
        f.write(str(sysfile))
        f.close()
        sys.exit(0)
    
    ################################################################################################
    # Output for CP2k
    if outputprogram == 'cp2k':
        docstring = StandardDocstring()
        cp2kinput = CP2KFile(cd, docstring)
        f = open(outputfile, outmode)
        f.write(str(cp2kinput))
        f.close()
        
    ################################################################################################
    # Output for CPMD
    if outputprogram == 'cpmd':
        docstring = StandardDocstring()
        cpmdinput = CPMDFile(cd, docstring)
        f = open(outputfile, outmode)
        f.write(str(cpmdinput))
        f.close()
        
    ################################################################################################
    # Output for Crystal09
    if outputprogram == 'crystal09':
        # Get documentation string
        docstring = StandardDocstring()
        # Get output file object
        crystal09file = Crystal09File(cd, docstring)
        crystal09file.spacegroupnr = cd.spacegroupnr
        crystal09file.a = cd.a
        crystal09file.b = cd.b
        crystal09file.c = cd.c
        crystal09file.alpha = cd.alpha
        crystal09file.beta  = cd.beta
        crystal09file.gamma = cd.gamma
        if cd.crystal_system() == "trigonal":
            # Get rhombohedral cell parameters
            if reducetoprim and abs(cd.gamma-120) < cd.coordepsilon:
                a = sqrt(3*cd.a**2 + cd.c**2)/3
                alpha = 360*asin(3/(2*sqrt(3+(cd.c/cd.a)**2)))/pi
                crystal09file.a = a
                crystal09file.b = a
                crystal09file.c = a
                crystal09file.alpha = alpha
                crystal09file.beta = alpha
                crystal09file.gamma = alpha
                crystal09file.trigonalsetting = "R"
        # Print to file
        f = open(outputfile, outmode)
        f.write(str(crystal09file))
        f.close()
        sys.exit(0)
    
    ################################################################################################
    # EMTO PROGRAMS
    if outputprogram == "emto" or outputprogram == "kgrn" or outputprogram == "kfcd" or outputprogram == "kstr" or outputprogram == "bmdl" or outputprogram == "shape":
        # Get job names
        if cd.HMSymbol == "":
            kstrjobnam = ref.cpd.replace(" ", "").replace("(","").replace(")","")
        else:
            kstrjobnam = cd.HMSymbol.replace("/","")
        kgrnjobnam = ref.cpd.replace(" ", "").replace("(","").replace(")","")
        # If the jobnames are long or blank, replace by the name of the cif file
        if len(kstrjobnam) > 30 and len(cif_file_name) < len(kstrjobnam)+4 or kstrjobnam == "":
            kstrjobnam = cif_file_name.replace(".cif","")
        if len(kgrnjobnam) > 30 and len(cif_file_name) < len(kgrnjobnam)+4 or kgrnjobnam == "":
            kgrnjobnam = cif_file_name.replace(".cif","")
        # Build docstring
        docstring = ref.compound+", "+ref.referencestring()
        # Document details of creation
        programdoc = "Generated by "+programname+" "+version+" "+datestring
        # Set the lattice number
        if cd.crystal_system() == "cubic":
            if cd.HMSymbol[0] == "F":
                latticenr = 2
            elif cd.HMSymbol[0] == "I":
                latticenr = 3
            else:
                latticenr = 1
        elif cd.crystal_system() == "hexagonal":
            latticenr = 4
        elif cd.crystal_system() == "tetragonal":
            if cd.HMSymbol == "I":
                latticenr = 6
            else:
                latticenr = 5
        elif cd.crystal_system() == "trigonal":
            latticenr = 7
        elif cd.crystal_system() == "orthorhombic":
            if cd.HMSymbol[0] == "A":
                latticenr = 9
            elif cd.HMSymbol[0] == "B":
                latticenr = 9
            elif cd.HMSymbol[0] == "C":
                latticenr = 9
            elif cd.HMSymbol[0] == "I":
                latticenr = 10
            elif cd.HMSymbol[0] == "F":
                latticenr = 11
            else:
                latticenr = 8
        elif cd.crystal_system() == "monoclinic":
            if cd.HMSymbol[0] == "A":
                latticenr = 13
            elif cd.HMSymbol[0] == "B":
                latticenr = 13
            elif cd.HMSymbol[0] == "C":
                latticenr = 13
            else:
                latticenr = 12
        else:
            # Triclinic and default
            latticenr = 14
    
    # Output for EMTO slope matrix program kstr
    if outputprogram == 'kstr' or outputprogram == 'emto':
        # Create directories
        try:
            os.mkdir('kstr')
        except OSError:
            pass
        try:
            os.mkdir('kstr/smx')
        except OSError:
            pass
        # Initialize BSTRFile
        kstrfile = KSTRFile(cd, docstring)
        kstrfile.jobnam = kstrjobnam
        if options.hardsphereradii:
            kstrfile.hardsphere = float(options.hardsphereradii)
        kstrfile.latticenr = latticenr
        kstrfile.a = cd.a
        kstrfile.b = cd.b
        kstrfile.c = cd.c
        kstrfile.alpha = cd.alpha
        kstrfile.beta = cd.beta
        kstrfile.gamma = cd.gamma
        # Document program
        kstrfile.programdoc = programdoc
        f = open("kstr/"+kstrjobnam+".dat","w")
        f.write(str(kstrfile))
        f.close()
        if outputprogram == "kstr":
            sys.exit(0)
    # Output for EMTO Madelung constant program bmdl
    if outputprogram == 'bmdl' or outputprogram == 'emto':
        # Create directories
        try:
            os.mkdir('bmdl')
        except OSError:
            pass
        try:
            os.mkdir('bmdl/mdl')
        except OSError:
            pass
        # Initialize BMDLFile
        bmdlfile = BMDLFile(cd, docstring)
        bmdlfile.jobnam = kstrjobnam
        bmdlfile.latticenr = latticenr
        bmdlfile.a = cd.a
        bmdlfile.b = cd.b
        bmdlfile.c = cd.c
        bmdlfile.alpha = cd.alpha
        bmdlfile.beta = cd.beta
        bmdlfile.gamma = cd.gamma
        # Document program
        bmdlfile.programdoc = programdoc
        f = open("bmdl/"+kstrjobnam+".dat","w")
        f.write(str(bmdlfile))
        f.close()
        if outputprogram == "bmdl":
            sys.exit(0)
    # Output for EMTO shape function program 'shape'
    if outputprogram == 'shape' or outputprogram == 'emto':
        # Create directories
        try:
            os.mkdir('shape')
        except OSError:
            pass
        try:
            os.mkdir('shape/shp')
        except OSError:
            pass
        # Initialize ShapeFile
        shapefile = ShapeFile(cd, docstring)
        shapefile.jobnam = kstrjobnam
        # Document program
        shapefile.programdoc = programdoc
        f = open("shape/"+kstrjobnam+".dat","w")
        f.write(str(shapefile))
        f.close()
        if outputprogram == "shape":
            sys.exit(0)
    # Output for EMTO main Greens function program 'kgrn'
    if outputprogram == 'kgrn' or outputprogram == 'emto':
        # Create directories
        try:
            os.mkdir('kgrn')
        except OSError:
            pass
        try:
            os.mkdir('kgrn/pot')
        except OSError:
            pass
        try:
            os.mkdir('kgrn/chd')
        except OSError:
            pass
        # Initialize KGRNFile
        kgrnfile = KGRNFile(cd, docstring)
        kgrnfile.jobnam = kgrnjobnam
        kgrnfile.kstrjobnam = kstrjobnam
        kgrnfile.latticenr = latticenr
        # Document program
        kgrnfile.programdoc = programdoc
        f = open("kgrn/"+kgrnjobnam+".dat","w")
        f.write(str(kgrnfile))
        f.close()
        if outputprogram == "kgrn":
            sys.exit(0)
    # Output for EMTO charge density program 'kfcd'
    if outputprogram == 'kfcd' or outputprogram == 'emto':
        # Create directories
        try:
            os.mkdir('kfcd')
        except OSError:
            pass
        # Initialize KFCDFile
        kfcdfile = KFCDFile(cd, docstring)
        kfcdfile.jobnam = kgrnjobnam
        kfcdfile.kstrjobnam = kstrjobnam
        # Document program
        kfcdfile.programdoc = programdoc
        f = open("kfcd/"+kgrnjobnam+".dat","w")
        f.write(str(kfcdfile))
        f.close()
        if outputprogram == "kfcd":
            sys.exit(0)        
    if outputprogram == "emto":
        sys.exit(0)
    
    ################################################################################################
    # Output for elk
    if outputprogram == 'elk':
        # Get documentation string
        docstring = StandardDocstring()
        # Get output file object
        geometryfile = ElkFile(cd, docstring)
        f = open(outputfile, outmode)
        f.write(str(geometryfile))
        f.close()
        sys.exit(0)
    
    # Output for exciting
    if outputprogram == 'exciting':
        # Get documentation string
        docstring = StandardDocstring()
        # Get output file object
        excitingfile = ExcitingFile(cd, docstring)
        f = open(outputfile, outmode)
        f.write(str(excitingfile))
        f.close()
        sys.exit(0)
    
    # Output for spacegroup.in for Elk/Exciting cell utility spacegroup
    if outputprogram == 'spacegroup':
        # Get documentation string
        docstring = StandardDocstring()
        # Get output file object
        spacegroupfile = SpacegroupFile(cd, docstring)
        angtobohr = 1e-10 * 4 * pi * 10973731.568527/7.2973525376e-3
        spacegroupfile.HermannMauguin = cd.HMSymbol.replace("/","")
        spacegroupfile.a = cd.a*angtobohr
        spacegroupfile.b = cd.b*angtobohr
        spacegroupfile.c = cd.c*angtobohr
        spacegroupfile.alpha = cd.alpha
        spacegroupfile.beta  = cd.beta
        spacegroupfile.gamma = cd.gamma
        if options.spacegroupsupercell:
            spacegroupfile.supercelldims = safe_matheval(options.spacegroupsupercell)
        # Print to file
        f = open(outputfile, outmode)
        f.write(str(spacegroupfile))
        f.close()
        #sys.exit(0)
    
    # Output for abacus originally from elk
    # !!!TODO: for general cifs  
    if outputprogram in ['abacus','abacus-lcao'] :
        print "    Name of CIF File:", os.path.split(options.file )[-1] 
        
        # Get documentation string
        docstring = StandardDocstring()
        # Get output file object
        #print "    kresolution",kresolution
        #geometryfile = AbacusFile(cd, docstring)
        geometryfile = AbacusFile(cd, docstring, kresolution=kresolution )

        geometryfile.cifname = os.path.split(options.file )[-1]
        #geometryfile.pspname = os.path.split(options.pspname )[-1]
        #print "    Name of CIF File:",geometryfile.cifname
        #print "    pspname:",geometryfile.pspname

        if options.deltav:
            geometryfile.deltav = safe_matheval(options.deltav)
        #if options.kptresol:
        #    geometryfile.kptresol = safe_matheval(options.kptresol)
        if options.pseudostringsuffix:
            geometryfile.pseudostringsuffix = options.pseudostringsuffix
            #print "    geometryfile.pseudostringsuffix = ", geometryfile.pseudostringsuffix
        if options.pspdir:
            geometryfile.pspdir = options.pspdir
            #print("safe_matheval->geometryfile.deltav", geometryfile.deltav)
        if options.energycutoff:
            geometryfile.energycutoff = safe_matheval(options.energycutoff)
        if outputprogram == 'abacus-lcao':
            geometryfile.uselcao = True
            #print "geometryfile.uselcao = ", geometryfile.uselcao
        if outputprogram == 'abacus':
            geometryfile.uselcao = False
        if options.orbfilepath:  
            geometryfile.orbfilepath = options.orbfilepath

        #print "\n1geometryfile.cell.lengthscale",geometryfile.cell.lengthscale
        #print   "1geometryfile.scale_newv",geometryfile.scale_newv, "\n"

        #print "outmode = ", outmode
        f = open("STRU", outmode)
        f.write(str(geometryfile))
        f.close() 
        #
        CellNatoms = geometryfile.cell.natoms()
        #CellNatoms = cd.numberOfAtoms 
        #CellNatoms = cd.natoms()
        CellScale  = geometryfile.scale_newv #inbohr
        CellVolume = geometryfile.cell.volume()*CellScale**3  #in bohr^3
        CellVolume = CellVolume * bohr_b_Ang**3  #in Ang^3

        print "    Cell Natoms = ", CellNatoms
        print "    Cell Scale  = ", CellScale  , " bohr "
        print "    Cell Volume = ", CellVolume , " Ang^3 "

        #print "2geometryfile.cell.lengthscale",geometryfile.cell.lengthscale
        #print "2geometryfile.scale_newv",geometryfile.scale_newv, "\n"


        # Output k-space mesh
        #nkpoints = int( (6750./geometryfile.cell.natoms() )**(1./3.) + 1)
        KPTa = geometryfile.kgrid[0]
        KPTb = geometryfile.kgrid[1]
        KPTc = geometryfile.kgrid[2]
        fKPT = open("KPT", 'w')
        fKPT.write(geometryfile.fKPTstring)
        fKPT.close() 
        KPTstring = ( "|+  KPTold: %7s"%geometryfile.cifname +" nkpts= %5d"%geometryfile.nkpoints_old 
                     +" ~%2d^3"%int((geometryfile.nkpoints_old)**(1./3)+1) 
                     +" kresol= %.3e /A"%(geometryfile.kresolution_old/bohr_b_Ang) 
                     +" kgrid_old~ [ %2d "%geometryfile.kgrid_old[0] 
                     +", %2d "%geometryfile.kgrid_old[1]
                     +", %2d ]"%geometryfile.kgrid_old[2] 
                     +" \n"
                     +"|+  KPTnow: %7s"%geometryfile.cifname +" nkpts= %5d"%geometryfile.nkpoints 
                     +" ~%2d^3"%int((geometryfile.nkpoints)**(1./3)+1) 
                     +" kresol= %.3e /A"%(geometryfile.kresolution/bohr_b_Ang) 
                     +" kgrid_now= [ %2d "%geometryfile.kgrid[0] 
                     +", %2d "%geometryfile.kgrid[1] 
                     +", %2d ]"%geometryfile.kgrid[2] 
                     ) 
        #print KPTstring

        # Output INPUT
        if outputfile :
            fINPUT = open(outputfile, 'w')
            fINPUT.write(geometryfile.fINPUTstring)
            fINPUT.close() 


        #sys.exit(0)
        
    ################################################################################################
    # Output for PWSCF (Quantum Espresso)
    if outputprogram == 'pwscf':
        docstring = StandardDocstring()
        pwscfinput = PWSCFFile(cd, docstring,kresolution=kresolution)
        #print(pwscfinput)
        if options.setupall:
            pwscfinput.setupall = True
            pwscfinput.kresolution = kresolution
        if options.pwscfpseudostring:
            pwscfinput.pseudostring = options.pwscfpseudostring
        if options.pwscfcart:
            pwscfinput.cartesian = True
        if options.pwscfcartvects:
            pwscfinput.cartesianlatvects = True
        if options.pwscfcartpos:
            pwscfinput.cartesianpositions = True
        if options.pwscfalatunits:
            pwscfinput.scaledcartesianpositions = True
        if options.pwscfatomicunits:
            pwscfinput.unit = "bohr"
        if options.deltav:
            pwscfinput.deltav = safe_matheval(options.deltav)
        if options.pseudostringsuffix:
            pwscfinput.pseudostringsuffix = options.pseudostringsuffix
        if options.pspdir:
            pwscfinput.pspdir = options.pspdir
            #print("safe_matheval->pwscfinput.deltav", pwscfinput.deltav)
        if options.energycutoff:
            pwscfinput.energycutoff = safe_matheval(options.energycutoff)
        #
        pwscfinput.cifname = os.path.split(options.file )[-1]
        pwscfinput.pspname = os.path.split(options.pspname )[-1]

        # 
        print "    pwscfinput.unit: ", pwscfinput.unit
        print "    cifname: ", pwscfinput.cifname
        print "    pspname: ", pwscfinput.pspname
        f = open(outputfile, outmode)
        f.write(str(pwscfinput))
        f.close()

        #calc CellScale (in bohr^3)
        if pwscfinput.unit == "bohr" :
            CellScale = pwscfinput.scale_newv  # in unit bohr 
        else:
            CellScale = pwscfinput.scale_newv / bohr_b_Ang  # in unit bohr 

        CellNatoms = pwscfinput.cell.natoms() 
        #CellNatoms = cd.numberOfAtoms 
        #CellNatoms = cd.natoms()

        #calc CellVolume (in Ang^3)
        CellVolume = cd.volume()*CellScale**3  # (in Ang^3) 
        if pwscfinput.unit == "bohr" :
            CellVolume = cd.volume() * pwscfinput.scale_newv**3 * bohr_b_Ang**3   # (in Ang^3) 
        else:
            CellVolume = cd.volume() * pwscfinput.scale_newv**3                   # (in Ang^3) 

        print "    CIF  Natoms = ", CellNatoms
        print "    CIF  Scale  = ", CellScale  , " bohr "
        print "    CIF  Volume = ", CellVolume , " Ang^3 "

        #CellVolume = CellVolume / pwscfinput.cell.natoms()  #(in A^3/atom)

        
    ################################################################################################
    # Output for Hutsepot
    if outputprogram == 'hutsepot':
        # Get documentation string
        docstring = StandardDocstring()
        # Get output file object
        hutsepotinput = HUTSEPOTFile(cd, docstring)
        f = open(outputfile,outmode)
        f.write(str(hutsepotinput))
        f.close()
    
    ################################################################################################
    # Output for FHI-AIMS
    if outputprogram == 'fhi-aims':
        # Get documentation string
        docstring = StandardDocstring()
        # Get output file object
        aimsinput = AIMSFile(cd, docstring)
        if options.aimscartesian:
            aimsinput.cartesian = True
        f = open(outputfile,outmode)
        f.write(str(aimsinput))
        f.close()
    
    ################################################################################################
    # Output for Fleur
    if outputprogram == 'fleur':
        # The first line with info from the CIF file and species order
        docstring = "Generated by "+programname+" "+version+" : "+ref.cpd+" ("+ref.compound+")"+" :  "+ref.referencestring()
        fleurinput = FleurFile(cd, docstring)
        f = open(outputfile, outmode)
        f.write(str(fleurinput))
        f.close()
        
    ################################################################################################
    # Output for mcsqs
    if outputprogram == 'mcsqs':
        # Get documentation string
        docstring = StandardDocstring()
        # Get output file object
        mcsqsinput = MCSQSFile(cd, docstring)
        f = open(outputfile,outmode)
        f.write(str(mcsqsinput))
        f.close()
    
    ################################################################################################
    # Output for MOPAC
    if outputprogram == 'mopac':
        docstring = StandardDocstring()
        mopacfirstline = ""
        mopacsecondline = ""
        mopacthirdline = ""
        mopacfreeze = -1
        mopacsetup = setupall
        if options.mopacfirstline:
            mopacsetup = True
            mopacfirstline = options.mopacfirstline
        if options.mopacsecondline:
            mopacsetup = True
            mopacsecondline = options.mopacsecondline
        if options.mopacthirdline:
            mopacsetup = True
            mopacthirdline = options.mopacthirdline
        if options.mopacfreeze:
            if options.mopacfreeze.upper() == 'F':
                mopacfreeze = 1
            elif options.mopacfreeze.upper() == 'T':
                mopacfreeze = 0
            else:
                sys.stderr.write("***Warning: I do not understand. --mopac-freeze-structure takes T (t) or F (f) as inputs.\n")
        mopacinput = MOPACFile(cd, docstring,setupall=mopacsetup,firstline=mopacfirstline,\
                               secondline=mopacsecondline,thirdline=mopacthirdline,freeze=mopacfreeze)
        f = open(outputfile, outmode)
        f.write(str(mopacinput))
        f.close()
        
    ################################################################################################
    # Output for TB-LMTO program ncol
    if outputprogram == 'ncol' or outputprogram == 'bstr':
        # Set up names for files. 
        if cd.HMSymbol == "":
            bstrjobnam = ref.cpd.replace(" ", "").replace("(","").replace(")","")
        else:
            bstrjobnam = cd.HMSymbol.replace("/","")
        jobnam = ref.cpd.replace(" ", "").replace("(","").replace(")","")
        # If the jobnames are long, replace by the name of the cif file
        if len(bstrjobnam) > 30 and len(cif_file_name) < len(bstrjobnam)+4:
            bstrjobnam = cif_file_name.replace(".cif","")
        if len(jobnam) > 10:
            if len(cif_file_name) < len(jobnam):
                jobnam = cif_file_name.replace(".cif","")
            else:
                jobnam = jobnam[0:9]
        # Build docstring
        docstring = ref.compound+", "+ref.referencestring()
        # Document details of creation
        programdoc = "Generated by "+programname+" "+version+" "+datestring
            
    # Output for TB-LMTO structure constant program bstr
    if outputprogram == 'bstr' or outputprogram == 'ncol':
        # Initialize BSTRFile
        bstrfile = BSTRFile(cd, docstring)
        bstrfile.jobnam = bstrjobnam
        bstrfile.a = cd.a
        bstrfile.b = cd.b
        bstrfile.c = cd.c
        # Document program
        bstrfile.programdoc = programdoc
        f = open(bstrjobnam+".dat","w")
        f.write(str(bstrfile))
        f.close()
        if outputprogram == 'bstr':
            sys.exit(0)
    
    if outputprogram == 'ncol':
        # Initialize ncolfile
        ncolfile = OldNCOLFile(cd, docstring)
        ncolfile.jobnam = jobnam
        ncolfile.bstrjobnam = bstrjobnam
        # Document program
        ncolfile.programdoc = programdoc
        f = open(jobnam+".dat","w")
        f.write(str(ncolfile))
        f.close()
        sys.exit(0)
    
        
    ################################################################################################
    # Output for VASP
    if outputprogram == 'vasp':
        # The first line with info from the CIF file and species order
        docstring =  "Generated by "+programname+" "+version
        if ref.database != "" and ref.databasecode != "":
            docstring += " from "+ref.databaseabbr[ref.database]+" reference: "+ref.databasecode
        docstring += ". "
        if len(ref.cpd) > len(ref.compound):
            docstring += ref.compound
        else:
            docstring += ref.cpd
        docstring += " :  "+ref.referencestring()+"."
        # Initialize the POSCARFile structure
        poscar = POSCARFile(cd, docstring, vca=options.vca)
        # Convert to cartesian coordinates if requested
        if options.vaspcartpos:
            poscar.printcartpos = True
        if options.vaspcartvecs:
            poscar.printcartvecs = True
        if options.vaspcart:
            poscar.printcartpos = True
            poscar.printcartvecs = True
        if options.vaspformat == "5":
            poscar.vasp5format = True
        if options.vaspselectivedyn:
            poscar.selectivedyn = True
        f = open(outputfile, outmode)
        f.write(str(poscar))
        f.close()
        # Print species order to screen if requested
        if options.vaspprintspcs:
            print poscar.SpeciesOrder()
        # Set up all files?
        if setupall:
            # POTCAR
            if options.vasppseudolib:
                lib = options.vasppseudolib
            else:
                lib = ""
            # Make selection of potcars
            if options.vasppppriority:
                prioritylist = options.vasppppriority.split(",")
                prioritylist.append("")
                print prioritylist
            else:
                try:
                    pl = os.environ['VASP_PP_PRIORITY']
                    prioritylist = pl.split(",")
                    prioritylist.append("")
                except:
                    prioritylist=["_d","_pv","_sv","","_h","_s"]
            if options.vaspencutfac:
                encutfac=float(options.vaspencutfac)
            else:
                encutfac=1.5
            potcarfile = POTCARFile(cd,directory=lib,vca=options.vca,prioritylist=prioritylist)
            f = open("POTCAR", "w")
            f.write(str(potcarfile))
            f.close()
            # KPOINTS
            docstring="Generated by cif2cell "+version+"."
            kpointsfile = KPOINTSFile(cd,docstring=docstring,kresolution=kresolution)
            f = open("KPOINTS", "w")
            f.write(str(kpointsfile))
            f.close()
            # INCAR
            incarfile = INCARFile(cd,docstring=docstring,vca=options.vca,prioritylist=prioritylist,encutfac=encutfac)
            f = open("INCAR", "w")
            f.write(str(incarfile))
            f.close()
        sys.exit(0)
        
    ################################################################################################
    # Output for RSPt
    if outputprogram == 'rspt':
        # Construct documentation string
        docstring = StandardDocstring()
        # Get file string and print to symt.inp/rspt.inp
        if options.newsymt or setupall:
            symtfile = SymtFile2(cd, docstring, kresolution=kresolution/angtobohr)
            # k-mesh generation etc
            if setupall:
                symtfile.setupall = True
            # spin polarization and relativity
            if options.rsptspinpol:
                symtfile.spinpol = True
            if options.rsptrelativistic:
                symtfile.relativistic = True
            if options.rsptmtradii:
                symtfile.mtradii = int(options.rsptmtradii)
            if options.rsptnospin:
                symtfile.forcenospin = True
        else:
            symtfile = SymtFile(cd, docstring)
        # spin axis 
        if options.rsptspinaxis:
            symtfile.spinaxis = Vector(safe_matheval(options.rsptspinaxis))
        if options.rsptpasswyckoff:
            symtfile.passwyckoff = True
        if options.rsptcartlatvects:
            symtfile.rsptcartlatvects = True
        #
        if options.exportlabels:
            symtfile.printlabels = True
        f = open(outputfile, outmode)
        f.write(str(symtfile))
        f.close()
        sys.exit(0)
        
    # Output for RSPt supercell generator cellgen
    if outputprogram == 'cellgen':
        # Construct documentation string
        docstring = StandardDocstring()
        # Initialize file object
        cellgenfile = CellgenFile(cd, docstring)
        if options.cellgenrefvec:
            cellgenfile.referencevector = safe_matheval(options.cellgenrefvec)
        if options.cellgensupercelldims:
            tmplist = safe_matheval(options.cellgensupercelldims)
            tmpmat = []
            for i in range(3):
                tmpmat.append([])
                for j in range(3):
                    tmpmat[i].append(0)
                tmpmat[i][i] = tmplist[i]
            cellgenfile.supercellmap = tmpmat
        if options.cellgenmap:
            tmpmat = safe_matheval(options.cellgenmap)
            cellgenfile.supercellmap = tmpmat
        # Print to file
        f = open(outputfile, outmode)
        f.write(str(cellgenfile))
        f.close()
        sys.exit(0)
    
    ################################################################################################
    # Output for Siesta
    if outputprogram == 'siesta':
        docstring = StandardDocstring()
        siestainput = SiestaFile(cd, docstring)
        f = open(outputfile, outmode)
        f.write(str(siestainput))
        f.close()
        
    ################################################################################################
    # Output for SPC file
    if outputprogram == 'spc':
        # The second comment line with info from the CIF file 
        docstring = StandardDocstring()
        # Initialize the SPCFile structure
        sysfile = SPCFile(cd, docstring)
        f = open(outputfile, outmode)
        f.write(str(sysfile))
        f.close()
        sys.exit(0)
    
    ################################################################################################
    # Output for xband
    if outputprogram == 'xband' or outputprogram == 'sprkkr':
        # The first line with info from the CIF file and species order
        docstring =  "Generated by "+programname+" "+version
        if ref.database != "" and ref.databasecode != "":
            docstring += " from "+ref.databaseabbr[ref.database]+" reference: "+ref.databasecode
        docstring += ". "
        if len(ref.cpd) > len(ref.compound):
            docstring += ref.compound
        else:
            docstring += ref.cpd
        docstring += " :  "+ref.referencestring()+"."
        # Initialize the file
        sysfile = XBandSysFile(cd, docstring)
        sysfile.filename = outputfile
        if options.sprkkrminangmom:
            sysfile.minangmom=int(options.sprkkrminangmom)
        f = open(outputfile, outmode)
        f.write(str(sysfile))
        f.close()
        sys.exit(0)
        
    ################################################################################################
    # Output for xyz file
    if outputprogram == 'xyz':
        # The second comment line with info from the CIF file 
        docstring =  "Generated by "+programname+" "+version
        if ref.database != "" and ref.databasecode != "":
            docstring += " from "+ref.databaseabbr[ref.database]+" reference: "+ref.databasecode
        docstring += ". "
        if len(ref.cpd) > len(ref.compound):
            docstring += ref.compound
        else:
            docstring += ref.cpd
        docstring += " :  "+ref.referencestring()+"."
        # Initialize the XYZFile structure
        if options.xyzatomicunits:
            cd.newunit("bohr")
        sysfile = XYZFile(cd, docstring)
        f = open(outputfile, outmode)
        f.write(str(sysfile))
        f.close()
        sys.exit(0)


if  __name__ == '__main__':
#    options.file='/home/nic/wszhang/eclipse_project/delta_dft/CIF_POT/CIF_elements/Al.cif'
#    options.program='abacus-lcao'
#    #options.program='pwscf'
#    options.pspdir = '/home/nic/wszhang/eclipse_project/delta_dft/CIF_POT/SG15_ONCV_PBE-1.0/'
#    options.pspname = 'Al_ONCV_PBE-1.0.upf' 
#    options.kresolution  = 0.0529177210903
#    options.energycutoff = '100'
#    #options.energycutoff = '80'
#    options.deltav = '1.0'
#    options.verbose = True

    #options.setupall = True 
    #options.pwscfatomicunits=True

    #options.pwscfpseudostring
    #options.abacusdeltav=1.1
    #print("\n")
    #print(options)
    #print(args)
    #print(options.append)
    #print(options.outputfile)

    print("    Running cif2cellmain")
    cif2cellmain() 


